home *** CD-ROM | disk | FTP | other *** search
/ Chip 1997 October / CD CHIP.ISO / Vari / Multimed / 95ITERAT / _SETUP.1 / Delta_z.cpp < prev    next >
Encoding:
C/C++ Source or Header  |  1997-01-13  |  27.2 KB  |  1,371 lines
  1. ///////////////////////////////////////////////
  2. // delta_z.cpp
  3.  
  4. #include "stdafx.h"
  5. #include "itriazon.h"
  6. #include "itriadoc.h"
  7. #include "itriavw.h"
  8. #include "external.h"
  9. #include "rw1.h"
  10. #include "rw2.h"
  11. #include "math.h"
  12. #include "convolut.h"
  13. #include "nthorder.h"
  14. #include "post.h"
  15. #include "orient.h"
  16. #include "filter12.h"
  17.  
  18. void CIterationsView::Delta_z(double dzx, double dzy)
  19. {
  20.     switch (nFilter)
  21.     {
  22.         case 1:    
  23.  
  24.             if (CRMAX > 2 || CIMAX > 2)
  25.             {
  26.                 if (CRMAX > CIMAX)
  27.                     dmax = CRMAX;
  28.                 else
  29.                     dmax = CIMAX;
  30.             }
  31.             else
  32.                 dmax = 2;
  33.  
  34.             delta5 = dmax/5;
  35.             
  36.             ntemp=(int)log(fabs(dzx*dzx+dzy*dzy));
  37.  
  38.             for (dlz=0 ; dlz<=dmax ; dlz+=delta5)
  39.             {    
  40.                 if (fabs(dzx) > dlz || fabs(dzy) > dlz) jrw+=ntemp;
  41.                 else                                                                        jrw-=ntemp;
  42.         }
  43.         break;
  44.     
  45.         case 2:    
  46.             if (CRMAX > 2 || CIMAX > 2)
  47.             {
  48.                 if (CRMAX > CIMAX)
  49.                     dmax = CRMAX;
  50.                 else
  51.                     dmax = CIMAX;
  52.             }
  53.             else
  54.                 dmax = 2;
  55.  
  56.             delta5 = dmax/5;
  57.  
  58.             ntemp = (int)(1/(sqrt(dzx*dzx+dzy*dzy)));
  59.             for (dlz=0 ; dlz<=dmax ; dlz+=delta5)
  60.             {    
  61.                 if (fabs(dzx) > dlz || fabs(dzy) > dlz) jrw-=ntemp;
  62.                 else                                                                        jrw+=ntemp;
  63.       }
  64.       break;
  65.       
  66.         case 3:    
  67.             if (CRMAX > 2 || CIMAX > 2)
  68.             {
  69.                 if (CRMAX > CIMAX)
  70.                     dmax = CRMAX;
  71.                 else
  72.                     dmax = CIMAX;
  73.             }
  74.             else
  75.                 dmax = 2;
  76.  
  77.             delta10 = dmax/10;
  78.  
  79.             for (dlz=0 ; dlz<=dmax ; dlz+=delta10)
  80.             {    
  81.                 if (fabs(dzx) > dlz || fabs(dzy) > dlz) jrw+=nFF;
  82.                 else                                                      jrw-=nFF;
  83.      }
  84.         break;
  85.         
  86.         case 4:  
  87.             for (dlz=0 ; dlz<=2 ; dlz+=.25)
  88.             {    
  89.                 if (fabs(dzx) > dlz && fabs(dzy) > dlz) jrw-=nFF;
  90.                 else                                                                         jrw+=nFF;
  91.         }
  92.         break;
  93.  
  94.         case 5:  
  95.             for (dlz=0 ; dlz<=2 ; dlz+=.25)
  96.             {    
  97.                 if (fabs(dzx) > dlz && fabs(dzy) > dlz) jrw-=nFF;
  98.                 else                                                                         jrw+=nFF;
  99.                 if (fabs(dzx) > dlz || fabs(dzy) > dlz) jrw-=nFF;
  100.                 else                                                                         jrw+=nFF;
  101.         }
  102.         break;
  103.  
  104.     case 6:      
  105.             if (dzx >    CRMID) jrw-=nFF;
  106.             else                         jrw+=nFF;
  107.             if (dzy    >    CIMID) jrw-=nFF;
  108.             else                         jrw+=nFF;
  109.           break;
  110.           
  111.         case 7:  
  112.  
  113.             //if (dzx*cy + cx*dzy >= 0)
  114.             //    jrw+=nFF;
  115.             //else
  116.             //    jrw-=nFF;
  117.  
  118.             if (dzx > CRMID && dzy < CIMID)     jrw+=nFF;
  119.                 else                                          jrw-=nFF;
  120.             if (dzx < CRMID && dzy > CIMID)     jrw+=nFF;
  121.                 else                                          jrw-=nFF;
  122.       break;
  123.             
  124.         case 8:  
  125.             if (dzx > CRMID && dzy > CIMID) jrw-=nFF;
  126.             else                                          jrw+=nFF;
  127.             if (dzx > CRMID && dzy > CIMID) jrw-=nFF;
  128.             else                                          jrw+=nFF;
  129.             if (dzx > CRMID || dzy > CIMID) jrw-=nFF;
  130.             else                                          jrw+=nFF;
  131.             if (dzx > CRMID || dzy > CIMID) jrw-=nFF;
  132.             else                                          jrw+=nFF;
  133.       break;
  134.  
  135.         case 9:  
  136.             
  137.             if (fabs(dzx) <= CRMID)               jrw+=nFF;
  138.             if (fabs(dzx) <= CRMID + .25)  jrw+=nFF;
  139.             if (fabs(dzx) <= CRMID + .50)  jrw+=nFF;
  140.             if (fabs(dzx) <= CRMID + .75)  jrw+=nFF;
  141.             if (fabs(dzx) <= CRMID + 1.0)  jrw+=nFF;
  142.             if (fabs(dzx) <= CRMID + 1.25) jrw+=nFF;
  143.             if (fabs(dzx) <= CRMID + 1.50) jrw+=nFF;
  144.             if (fabs(dzx) <= CRMID + 1.75) jrw+=nFF;
  145.             if (fabs(dzx) <= CRMID + 2.00) jrw+=nFF;
  146.  
  147.             if (fabs(dzy) <= CIMID)                 jrw+=nFF;
  148.             if (fabs(dzy) <= CIMID + .25)  jrw+=nFF;
  149.             if (fabs(dzy) <= CIMID + .50)  jrw+=nFF;
  150.             if (fabs(dzy) <= CIMID + .75)  jrw+=nFF;
  151.             if (fabs(dzy) <= CIMID + 1.0)  jrw+=nFF;
  152.             if (fabs(dzy) <= CIMID + 1.25) jrw+=nFF;
  153.             if (fabs(dzy) <= CIMID + 1.50) jrw+=nFF;
  154.             if (fabs(dzy) <= CIMID + 1.75) jrw+=nFF;
  155.             if (fabs(dzy) <= CIMID + 2.00) jrw+=nFF;
  156.  
  157.             break;
  158.  
  159.         case 10:  
  160.             dmax = fabs(CRMAX - CRMIN);
  161.  
  162.             if (fabs(dzx) <=     CRMID +      0) jrw+=nFF;
  163.             if (fabs(dzx) <=  CRMID + dmax/2) jrw+=nFF;
  164.             if (fabs(dzx) <=  CRMID + dmax/3) jrw+=nFF;
  165.             if (fabs(dzx) <=  CRMID + dmax/4) jrw+=nFF;
  166.             if (fabs(dzx) <=  CRMID + dmax/5) jrw+=nFF;
  167.             if (fabs(dzx) <=  CRMID + dmax/6) jrw+=nFF;
  168.             if (fabs(dzx) <=  CRMID + dmax/7) jrw+=nFF;
  169.             if (fabs(dzx) <=  CRMID + dmax/8) jrw+=nFF;
  170.             if (fabs(dzx) <=  CRMID + dmax/9) jrw+=nFF;
  171.  
  172.             dmax = fabs(CIMAX - CIMIN);
  173.  
  174.             if (fabs(dzy) <=     CIMID +      0) jrw+=nFF;
  175.             if (fabs(dzy) <=  CIMID + dmax/2) jrw+=nFF;
  176.             if (fabs(dzy) <=  CIMID + dmax/3) jrw+=nFF;
  177.             if (fabs(dzy) <=  CIMID + dmax/4) jrw+=nFF;
  178.             if (fabs(dzy) <=  CIMID + dmax/5) jrw+=nFF;
  179.             if (fabs(dzy) <=  CIMID + dmax/6) jrw+=nFF;
  180.             if (fabs(dzy) <=  CIMID + dmax/7) jrw+=nFF;
  181.             if (fabs(dzy) <=  CIMID + dmax/8) jrw+=nFF;
  182.             if (fabs(dzy) <=  CIMID + dmax/9) jrw+=nFF;
  183.             break;
  184.  
  185.         case 11:  
  186.             if (CRMAX > 2 || CIMAX > 2)
  187.             {
  188.                 if (CRMAX > CIMAX)
  189.                     dmax = CRMAX;
  190.                 else
  191.                     dmax = CIMAX;
  192.             }
  193.             else
  194.                 dmax = 2;
  195.  
  196.             delta10 = dmax/10;
  197.  
  198.             for (dlz=0 ; dlz<=dmax ; dlz+=delta10)
  199.             {    
  200.                 if (fabs(dzx) <=  dlz) jrw+=nFF;
  201.                 if (fabs(dzy) <=  dlz) jrw+=nFF;
  202.             }
  203.             break;
  204.  
  205.         case 12:  // Initialize c = (c^n).clog;
  206.             break;
  207.  
  208.         case 15:  //
  209.             if (dzx_save < fabs(dzx))
  210.             {
  211.                 dzx_save = fabs(dzx);
  212.                 jrw+=nFF;
  213.             }
  214.             if (dzy_save < fabs(dzy))
  215.             {
  216.                 dzy_save = fabs(dzy);
  217.                 jrw+=nFF;
  218.             }
  219.             break;
  220.  
  221.         case 16:  // 
  222.             if (dzx_save == 0)
  223.                 dzx_save = 99;
  224.             if (dzy_save == 0)
  225.                 dzy_save = 99;
  226.  
  227.             if (dzx_save > fabs(dzx))
  228.             {
  229.                 dzx_save = fabs(dzx);
  230.                 jrw+=nFF;
  231.             }
  232.             if (dzy_save > fabs(dzy))
  233.             {
  234.                 dzy_save = fabs(dzy);
  235.                 jrw+=nFF;
  236.             }
  237.             break;
  238.  
  239.         case 17:  //
  240.             if (dzx_save == 0)
  241.             {
  242.                 bPositiveX=FALSE;
  243.                 bPositiveY=FALSE;
  244.             }
  245.  
  246.             if (fabs(dzx) >= dzx_save)
  247.             {
  248.                 // Count the Positive Slope Change
  249.                 if (!bPositiveX)
  250.                 {
  251.                     jrw+=nFF;
  252.                     bPositiveX = TRUE;
  253.                 }
  254.             }
  255.             else
  256.             {
  257.                 // Count the Negative Slope Change
  258.                 if (bPositiveX)
  259.                 {
  260.                     jrw+=nFF;
  261.                     bPositiveX = FALSE;
  262.                 }
  263.             }            
  264.  
  265.             if (fabs(dzy) >= dzy_save)
  266.             {
  267.                 // Count the Positive Slope Change
  268.                 if (!bPositiveY)
  269.                 {
  270.                     jrw+=nFF;
  271.                     bPositiveY=TRUE;
  272.                 }
  273.             }
  274.             else
  275.             {
  276.                 // Count The Negative Slope Change
  277.                 if (bPositiveY)
  278.                 {
  279.                     jrw+=nFF;
  280.                     bPositiveY=FALSE;
  281.                 }
  282.             }
  283.             
  284.             dzx_save = fabs(dzx);
  285.             dzy_save = fabs(dzy);
  286.  
  287.             break;
  288.  
  289.         case 18:  //
  290.             if (dzx_save == 0)
  291.             {
  292.                 bPositiveX=FALSE;
  293.                 bPositiveY=FALSE;
  294.                 z1 = z;
  295.                 dzx_save = 99;
  296.                 jrw = NMAX/4;
  297.                 return;
  298.  
  299.             }
  300.  
  301.             //if (z.real() >= z1.real() || z.imag() >= z1.imag())
  302.             if (fabs(z.real()) >= fabs(z1.real()) || fabs(z.imag()) >= fabs(z1.imag()))
  303.             //if (z.squares() >= z1.squares())
  304.             {
  305.                 // Count the Positive Slope No Change
  306.                 if (!bPositiveX)
  307.                 {
  308.                     bPositiveX = TRUE;
  309.                     if (jrw-nFF > 0)
  310.                         jrw-=nFF;
  311.                 }
  312.                 else
  313.                 {
  314.                     if (jrw+(nFF*2) < NMAX)
  315.                         jrw+=(nFF*2);
  316.                 }
  317.             }
  318.             else
  319.             {
  320.                 // Count the Negative Slope No Change
  321.                 if (bPositiveX)
  322.                 {
  323.                     bPositiveX = FALSE;
  324.                     if (jrw-nFF > 0)
  325.                         jrw-=nFF;
  326.                 }
  327.                 else
  328.                 {
  329.                     if (jrw+(nFF*2) < NMAX)
  330.                         jrw+=(nFF*2);
  331.                 }
  332.             }            
  333.  
  334.             z1 = z;
  335.  
  336.             break;
  337.  
  338.         case 19:  
  339.             dzx_save += fabs(dzx);
  340.             dzy_save += fabs(dzy);
  341.             
  342.             break;
  343.  
  344.         case 20:  
  345.             if (dzx < .01 && dzx > -.01)
  346.             {
  347.                 if (dzx < 0)
  348.                     dzx = -.01;
  349.                 else
  350.                     dzx =  .01;
  351.             }
  352.  
  353.             if (dzy < .01 && dzy > -.01)
  354.             {
  355.                 if (dzy < 0)
  356.                     dzy = -.01;
  357.                 else
  358.                     dzy = .01;
  359.             }
  360.  
  361.             dzx_save += log(fabs(dzx));
  362.             dzy_save += log(fabs(dzy));
  363.  
  364.             break;
  365.  
  366.         case 23:
  367.  
  368.             k = i;
  369.             r = dzx*dzx+dzy*dzy;
  370.             if (k == 0) { k1 = 1; rz = r; }
  371.             else if (k == 1)
  372.             {
  373.                 if (r > rz) rzflag = 1;
  374.                 else rzflag = -1;
  375.                 rz = r;
  376.             }
  377.             else
  378.             {
  379.                 if (r > rz && rzflag != 1)
  380.                 {
  381.                     jrw+=nFF;
  382.                     rzflag = 1;
  383.                 }
  384.                 else 
  385.                 if (r < rz && rzflag != -1)
  386.                 {
  387.                     jrw+=nFF;
  388.                     rzflag = -1;
  389.                 }
  390.                 rz = r;
  391.             }
  392.             break;
  393.         
  394.         case 24:
  395.             z.set_real(z.real()-dF/sin(z.real()));
  396.             z.set_imag(z.imag()-dF/sin(z.imag()));
  397.  
  398.             break;
  399.         
  400.         case 25:
  401.             if (z.real() >= 0 && z.imag() >= 0) // Top,Right
  402.                 z = z-cmplx(dF,dF)/z.csin();
  403.             else
  404.             if (z.real() >= 0 && z.imag() < 0)  // Bottom, Right
  405.                 z = z-cmplx(dF,-dF)/z.csin();
  406.             else
  407.             if (z.real() < 0 && z.imag() < 0)        // Bottom, Left
  408.                 z = z-cmplx(-dF,-dF)/z.csin();
  409.             else
  410.                 z = z-cmplx(-dF,+dF)/z.csin();    // Top, Left
  411.  
  412.             break;
  413.                 
  414.         case 30:
  415.             // Counting
  416.  
  417.             if (fabs(dzx)>fabs(dzy))
  418.                     cNMAX.set_real(cNMAX.real() - fabs(dzx) - fabs(dzx_save));
  419.                 else
  420.                     cNMAX.set_imag(cNMAX.imag() - fabs(dzy) - fabs(dzy_save));
  421.  
  422.             dzx_save = dzx;
  423.             dzy_save = dzy;
  424.  
  425. //                cNMAX = z+cmplx((c.real())*(CRMAX-CIMIN), (c.imag())*(CIMAX-CIMIN));
  426.  
  427. //                cNMAX += cmplx(dzx+c.real()/(CRMAX-CIMIN), dzy+c.imag()/(CIMAX-CIMIN));
  428.  
  429. //            cNMAX = cNMAX+cmplx(1/dzx+c.real(),1/dzy+c.imag());
  430.             
  431. //            cNMAX = cNMAX+cmplx(1/dzx,1/dzy)+c;
  432.  
  433. //            cNMAX = cNMAX+cmplx(1/dzx,1/dzy)+c;
  434.             break;        
  435.         
  436.  
  437.         ////////////////////////////
  438.         // Category 2 filters
  439.         case 130:  //         
  440.             //if (fabs(dzx) < dStrands || fabs(dzy) < dStrands)
  441.                 //i3++;
  442.  
  443.             //rr = (sin(fabs(1 + cx_hi - cx_lo)) * sin(fabs(1 + cy_hi - cy_lo))) * 40;
  444.  
  445.             if (i <= 1)
  446.             {
  447.                 x_rmin = dzx;
  448.                 x_rmax = dzx;
  449.                 y_rmin = dzy;
  450.                 y_rmax = dzy;
  451.             }    
  452.             else
  453.             {    
  454.                 if (x_rmin > dzx)
  455.                     x_rmin = dzx;
  456.                 if (x_rmax < dzx)
  457.                     x_rmax = dzx;
  458.  
  459.                 if (y_rmin > dzy)
  460.                     y_rmin = dzy;
  461.                 if (y_rmax < dzy)
  462.                     y_rmax = dzy;
  463.             }
  464.  
  465.             break;        
  466.  
  467.         case 131:
  468.             if ((fabs(dzx) <= dStrands_HI && fabs(dzx) > dStrands_LO) || (fabs(dzy) <= dStrands_HI && fabs(dzy) > dStrands_LO))        
  469.              rr+=log(dzx*dzx+dzy*dzy)*13;  // Don't mess with this
  470.  
  471.             if (fabs(dzx) < dStrands || fabs(dzy) < dStrands)
  472.                 i3++;
  473.  
  474.             break;
  475.  
  476.         case 132:
  477.             if (dzx*dzx+dzy*dzy < limit)
  478.                 i3++;
  479.             break;
  480.  
  481.         case 133:
  482.             if (dzx < dStrands && dzx > -dStrands)
  483.             {
  484.                 if (dzx < 0)
  485.                     dzx = -dStrands;
  486.                 else
  487.                     dzx =  dStrands;
  488.             }
  489.  
  490.             if (dzy < dStrands && dzy > -dStrands)
  491.             {
  492.                 if (dzy < 0)
  493.                     dzy = -dStrands;
  494.                 else
  495.                     dzy = dStrands;
  496.             }
  497.  
  498.             dzx_save += log(fabs(dzx));
  499.             dzy_save += log(fabs(dzy));
  500.             
  501.             if ((fabs(dzx) <= dStrands_HI && fabs(dzx) > dStrands_LO) || (fabs(dzy) <= dStrands_HI && fabs(dzy) > dStrands_LO))        
  502.              rr+=log(dzx*dzx+dzy*dzy)*13;  // Don't mess with this            
  503.             
  504.             if (dzx*dzx+dzy*dzy < limit)
  505.                 i3++;
  506.  
  507.             break;
  508.  
  509.         case 134:
  510.            rr+=log(dzx*dzx+dzy*dzy)*10;        
  511.             if (fabs(dzx) > dStrands || fabs(dzy) > dStrands)
  512.                 i+=nFF;
  513.             //i%=NMAX;
  514.             break;
  515.             
  516.         case 135:
  517.             if (fabs(dzx) < limit || fabs(dzy) < limit)
  518.              rr+=log(dzx*dzx+dzy*dzy)*10;
  519.             else
  520.                 i+=nFF;
  521.             //i%=NMAX;
  522.             break;
  523.  
  524.         case 136:
  525.             rr+=log(dzx*dzx+dzy*dzy)*10;        
  526.             if (fabs(dzx) > limit || fabs(dzy) > limit)
  527.                 i+=nFF;
  528.             //i%=NMAX;
  529.             break;
  530.  
  531.         case 137:
  532.             if (fabs(dzx) < limit || fabs(dzy) < limit)
  533.              rr+=log(dzx*dzx+dzy*dzy)*10;
  534.             else
  535.                 i+=nFF;
  536.             //i%=NMAX;
  537.             break;
  538.  
  539.         case 138:
  540.             rr+=log(dzx*dzx+dzy*dzy)*10;        
  541.             if (fabs(1/dzx*dzy) > limit)
  542.                 i+=nFF;
  543.             //i%=NMAX;
  544.             break;
  545.  
  546.         case 139:
  547.             rr+=log(dzx*dzx+dzy*dzy)*10;        
  548.             if (fabs(dzx) > limit || fabs(dzy) > limit)
  549.                 i+=nFF;
  550.             //i%=NMAX;
  551.  
  552.             dzx_save = dzx;
  553.             dzy_save = dzy;
  554.             break;
  555.  
  556.         case 140:
  557.             rr+=log(dzx*dzx+dzy*dzy)*2;        
  558.             if (fabs(dzx) > limit || fabs(dzy) > limit)
  559.                 i+=nFF;
  560.             //i%=NMAX;
  561.             
  562.             dzx_save = dzx;
  563.             dzy_save = dzy;
  564.             break;
  565.  
  566.         case 141:
  567.             rr+=log(dzx*dzx+dzy*dzy)*nBay100;        
  568.             //if (fabs(dzx) > limit || fabs(dzy) > limit)
  569.                 //i+=nFF;
  570.             //i%=NMAX;
  571.             
  572.             dzx_save = dzx;
  573.             dzy_save = dzy;
  574.             
  575.             break;
  576.  
  577.         case 142:
  578.             rr+=log(dzx*dzx+dzy*dzy);        
  579.             if (fabs(dzx) > limit || fabs(dzy) > limit)
  580.                 i+=nFF;
  581.             //i%=NMAX;
  582.                         
  583.             dzx_save = dzx;
  584.             dzy_save = dzy;
  585.             break;
  586.  
  587.         case 143:
  588.             rr+=log(dzx*dzx+dzy*dzy)*8;        
  589.             if (fabs(dzx) > limit || fabs(dzy) > limit)
  590.                 i+=nFF;
  591.             //i%=NMAX;
  592.                         
  593.             dzx_save = dzx;
  594.             dzy_save = dzy;
  595.             break;
  596.  
  597.         case 144:
  598.             if (fabs(dzx) < dStrands || fabs(dzy) < dStrands)
  599.              //rr+=log(dzx*dzx+dzy*dzy)*3;  // Don't mess with this
  600.                 rr+=nFF;
  601.  
  602.             dzx_save = dzx;
  603.             dzy_save = dzy;
  604.             break;
  605.  
  606.         case 145:
  607.             if (fabs(dzx) < dStrands || fabs(dzy) < dStrands)
  608.              //rr+=log(dzx*dzx+dzy*dzy)*3;  // Don't mess with this
  609.                 rr+=nFF;
  610.  
  611.             dzx_save = dzx;
  612.             dzy_save = dzy;
  613.             break;
  614.         
  615.         case 146:
  616.             rr+=log(dzx*dzx+dzy*dzy)*(1+nBay100);        
  617.             
  618.             //if (fabs(dzx) > limit || fabs(dzy) > limit)
  619.                 //i+=nFF;
  620.             
  621.             //dzx_save = dzx;
  622.             //dzy_save = dzy;
  623.             
  624.             break;
  625.  
  626.         case 147:
  627.             rr+=log(dzx*dzx+dzy*dzy)*(1+nBay100);        
  628.             
  629.             //if (fabs(dzx) > limit || fabs(dzy) > limit)
  630.                 //i+=nFF;
  631.             
  632.             dzx_save = dzx;
  633.             dzy_save = dzy;
  634.             
  635.             break;
  636.                 
  637.         default:
  638.             break;
  639.     }
  640.     
  641.     if (bGeometry)
  642.     {
  643.         if (bFilter14)    // Swap real and imaginary components
  644.         {
  645.             z.set_real(dzy);
  646.             z.set_imag(dzx);
  647.         }        
  648.                 
  649.         if (bFilter21)        // Double
  650.         {
  651.             if (dzx >= 0)
  652.                 z.set_real(z.real()-2);
  653.             else
  654.                 z.set_real(z.real()+2);
  655.         }                
  656.  
  657.         if (bFilter22)        // Quad
  658.         {
  659.             if (dzx >= 0)
  660.              z.set_real(z.real()-2);
  661.             else
  662.                 z.set_real(z.real()+2);
  663.                     
  664.             if (dzy >= 0)
  665.                 z.set_imag(z.imag()-2);
  666.             else
  667.                 z.set_imag(z.imag()+2);
  668.         }                
  669.         
  670.         if (bFilter26)  // z^n
  671.             z=z^cFilter26;
  672.  
  673.         if (bFilter31) 
  674.             z=z*z.csin();
  675.  
  676.         if (bFilter32) 
  677.             z=z*z.ccos();
  678.     
  679.         if (bFilter33) 
  680.             z=z*z.clog();
  681.  
  682.         if (bFilter34) 
  683.             z=z*tangent(z);
  684.  
  685.         if (bFilter35) 
  686.             z=z*sinh(z);
  687.  
  688.         if (bFilter36) 
  689.             z=z*asin(z);
  690.     
  691.         if (bFilter37) 
  692.             z=z*acos(z);
  693.  
  694.         if (bFilter38) 
  695.             z=z*arctan(z);
  696.     }
  697.  
  698.     if (bMFilter)
  699.     {
  700.         // Make sure i is within limits
  701.         if (i >= NMAX || i < 0)
  702.             i = NMAX - 1;
  703.  
  704.         // Store the data into the temp arrays for Standard Deviation
  705.         // and Fractal Dimension calculations
  706.         if (nPreFractal)
  707.         {
  708.             pXTemp[i] = z.real();
  709.             pYTemp[i] = z.imag();    
  710.         }
  711.         else
  712.         {
  713.             pXTemp[i] = dzx;
  714.             pYTemp[i] = dzy;    
  715.         }
  716.     }
  717.  
  718.     if (bOrbits)
  719.     {
  720.         if (px == rb_center_x && py == rb_center_y)
  721.         {
  722.             CIterationsDoc* pDoc = GetDocument();
  723.  
  724.             pDoc->lpCt[i] = i;
  725.             pDoc->lpX1[i] = dzx;
  726.             pDoc->lpY1[i] = dzy;
  727.             
  728.             //return;
  729.  
  730.         }
  731.     }
  732.  
  733.     if (bQuickMode && nFilter != 29 )
  734.  
  735.     {
  736.         if (nDistortion <= 69 && nDistortion >= 61)
  737.             // Don't apply to Newton Type Fractals 61 - 69
  738.             return;
  739.  
  740.         if (fabs(dzx-dzx_save_quick) < dF28 && fabs(dzy-dzy_save_quick) < dF28)
  741.         {
  742.             if (bOrbits)
  743.             {
  744.                 for (i2 = i ; i2<NMAX ; i2++)
  745.                 {
  746.                     // fill remaining
  747.                     //pXTemp[i2] = dzx;
  748.                     //pYTemp[i2] = dzy;
  749.                     
  750.                     CIterationsDoc* pDoc = GetDocument();
  751.  
  752.                     pDoc->lpCt[i2] = i2;
  753.                     pDoc->lpX1[i2] = dzx;
  754.                     pDoc->lpY1[i2] = dzy;
  755.  
  756.                     //AfxMessageBox("Q_1");
  757.                 }
  758.             }
  759.  
  760.             bQuickAbort = TRUE;
  761.             iQuickAbort = i;
  762.  
  763.             i = NMAX;
  764.         }            
  765.         else
  766.             bQuickAbort = FALSE;
  767.  
  768.         //if (bQuickAbort && bOrbits)
  769.         //    AfxMessageBox("Q_2");
  770.         
  771.         dzx_save_quick = dzx;
  772.         dzy_save_quick = dzy;
  773.     }
  774. }    
  775.  
  776. void CIterationsView::Filter_Complete()
  777. {
  778.     if (i >= NMAX-1)
  779.         b_MAX = TRUE;
  780.     else
  781.         b_MAX = FALSE;
  782.     
  783.     if (bQuickMode && nFilter != 29)
  784.     {
  785.         if (bQuickAbort)
  786.             i = NMAX - 1;
  787.             //            i = iQuickAbort;
  788.  
  789.     }
  790.         
  791.     if (bOrbits)
  792.     {
  793.         if (px == rb_center_x && py == rb_center_y)
  794.         {
  795.             CIterationsDoc* pDoc = GetDocument();
  796.             pDoc->i_sav4orb = i;
  797.  
  798.             if (i < 0 || i > NMAX)
  799.             {
  800.                 char cstr[81];
  801.                 wsprintf(cstr,"Cannot Plot i=%d, Min=2, Max=%d",i, NMAX-1);
  802.                 AfxMessageBox(cstr);            
  803.                 return;
  804.             }
  805.  
  806. //            if (i < NMAX - 2)
  807. //            {
  808. //                i++;
  809. //                pDoc->i_sav4orb = i;
  810. //
  811. //                // Get Last bailout values
  812. //                pDoc->lpCt[i] = i;
  813. //                pDoc->lpX1[i2] = z.real();
  814. //                pDoc->lpY1[i2] = z.imag();
  815. //
  816. //            }
  817.  
  818.             //pDoc->lpCt[0] = 0;
  819.             //pDoc->lpX1[0] = 0;
  820.             //pDoc->lpY1[0] = 0;
  821.                         
  822.             // Initialize the range    with the first value
  823.             x_rmin = pDoc->lpX1[0];    // x range min
  824.             x_rmax = pDoc->lpX1[0];    // x range max
  825.             y_rmin = pDoc->lpY1[0];    // y range min
  826.             y_rmax = pDoc->lpY1[0];    // y range max
  827.  
  828.             for (i2 = 0 ; i2 <= i ; i2++)
  829.             {
  830.                 
  831.                 // Create range array
  832.                 // pDoc->lpR1[i2] = pDoc->lpX1[i2] + pDoc->lpY1[i2];
  833.                 if (pDoc->lpX1[i2] > x_rmax)
  834.                     x_rmax = pDoc->lpX1[i2];
  835.                 if (pDoc->lpX1[i2] < x_rmin)
  836.                     x_rmin = pDoc->lpX1[i2];
  837.  
  838.                 if (pDoc->lpY1[i2] > y_rmax)
  839.                     y_rmax = pDoc->lpY1[i2];
  840.                 if (pDoc->lpY1[i2] < y_rmin)
  841.                     y_rmin = pDoc->lpY1[i2];
  842.  
  843.                 if (i2 % 2 == 0)
  844.                 {
  845.                     // Even max
  846.                     if (x_rmax > y_rmax)
  847.                         pDoc->lpR1[i2] = x_rmax;
  848.                     else
  849.                         pDoc->lpR1[i2] = y_rmax;
  850.                 }
  851.                 else
  852.                 {
  853.                     // Odd min
  854.                     if (x_rmin < y_rmin)
  855.                         pDoc->lpR1[i2] = x_rmin;
  856.                     else
  857.                         pDoc->lpR1[i2] = y_rmin;
  858.                 }
  859.             }
  860.  
  861. //            if (x_rmin < y_rmin)
  862. //                pDoc->lpR1[i2-2] = x_rmin;
  863. //            else
  864. //                pDoc->lpR1[i2-2] = y_rmin;
  865. //
  866. //            if (x_rmax > y_rmax)
  867. //                pDoc->lpR1[i2-1] = x_rmax;
  868. //            else
  869. //                pDoc->lpR1[i2-1] = y_rmax;
  870.  
  871.             if (x_rmin < y_rmin)
  872.                 pDoc->lpR1[0] = x_rmin;
  873.             else
  874.                 pDoc->lpR1[0] = y_rmin;
  875.  
  876.             if (x_rmax > y_rmax)
  877.                 pDoc->lpR1[1] = x_rmax;
  878.             else
  879.                 pDoc->lpR1[1] = y_rmax;
  880.             
  881.             pDoc->UpdateAllGraphs();
  882.         }
  883.     }    
  884.  
  885.  
  886.     switch(nFilter)
  887.     {
  888.         case 14:
  889.             return;
  890.     
  891.         case 15:
  892.         case 16:
  893.         case 17:
  894.     
  895.             i = jrw;
  896.             i = abs(i);
  897.             i %= NMAX;
  898.             return;
  899.  
  900.         case 18:
  901.  
  902.             // An Inside Filter
  903.             if (i >= NMAX-1)
  904.                 i = jrw;
  905.             else
  906.                 i += jrw;
  907.     
  908.             i = abs(i);
  909.             i %= NMAX;
  910.             return;
  911.  
  912.         case 19:
  913.             if (i != 0) 
  914.                 i = (int) (((dzx_save + dzy_save)/(double)(i*2)) * 10 * nBay100);
  915.             i = abs(i);
  916.             if (i >= NMAX) i %= (NMAX-1);
  917.             return;    
  918.  
  919.         case 20:
  920.             if (dzx_save >= dzy_save)
  921.                 i = (int) ((dzx_save - dzy_save) * nBay100);
  922.             else
  923.                 i = (int) ((dzy_save - dzx_save) * nBay100);
  924.         
  925.             i = abs(i);
  926.             if (i >= NMAX) i %= NMAX;
  927.             return;    
  928.  
  929.         case 23:
  930.             // An Inside Filter
  931.             if (i >= NMAX-1)
  932.                 i = jrw;
  933.             else
  934.                 i += jrw;
  935.  
  936.             i = abs(i);
  937.             i %= NMAX;
  938.             return;
  939.  
  940.         case 30:
  941.             i = (int) fabs((cNMAX.squares() * nBay100));
  942.             if (i >= NMAX) i %= (NMAX-1);
  943.             return;    
  944.  
  945.         case 130:
  946.             //i = i3;
  947.             //i = (int) ((sin(fabs(1 + x_rmax - x_rmin)) * sin(fabs(1 + y_rmax - y_rmin))) * 40);
  948.             i = (int) (atan(fabs(x_rmax * x_rmin)/fabs(y_rmax * y_rmin)) * 40);
  949.  
  950.             return;            
  951.  
  952.         case 131:
  953.             if (i >= NMAX-1)
  954.                 i = i3;
  955.             else
  956.                 i = (int)rr;
  957.             return;
  958.  
  959.         case 132:
  960.             i = i3;
  961.             return;
  962.  
  963.         case 133:
  964.             if (rr != 0 && i3 != 0)
  965.                 i = (int)(fabs(dzy_save - dzx_save) + fabs(rr + i3));
  966.             else
  967.                 if (rr != 0)
  968.                     i = (int)(fabs(dzy_save - dzx_save) + fabs(rr));
  969.             else
  970.                 if (i3 != 0)
  971.                     i = (int)fabs(dzy_save - dzx_save) + i3;
  972.             else
  973.                 i = (int)fabs(dzy_save - dzx_save);
  974.             return;
  975.  
  976.         case 134:
  977.             rr = (double) ((int)rr%(NMAX/2));
  978.             if ((int)(z.real()*2)%2 == 0 || (int)(z.imag()*2)%2 == 0)
  979.                 i = (int)rr;
  980.             else
  981.                 i = (int)rr + (NMAX/2);
  982.             return;
  983.  
  984.         case 135:
  985.             rr = (double) ((int)rr%(NMAX/2));
  986.             if ((int)(z.real()*2)%2 == 0 || (int)(z.imag()*2)%2 == 0)
  987.                 i = (int)rr;
  988.             else
  989.                 i = (int)rr + (NMAX/2);
  990.             return;
  991.  
  992.         case 136:
  993.             i = (int)rr;
  994.             return;
  995.  
  996.         case 137:
  997.             i = (int)rr;
  998.             return;
  999.  
  1000.         case 138:
  1001.             i = (int)rr;
  1002.             return;
  1003.  
  1004.         case 139:
  1005.             i = (int)(atan(dzx_save/dzy_save) * 80); 
  1006.             return;
  1007.  
  1008.         case 140:
  1009.             i = (int)(rr + atan(fabs(dzx_save/dzy_save)) * 30); 
  1010.             return;
  1011.  
  1012.         case 141:
  1013.             //i += (int)(rr); 
  1014.             i = (int)rr;
  1015.             
  1016.             //i = (int)(rr + atan(fabs(dzx_save/dzy_save)) * 30); 
  1017.  
  1018. //            i = (int)(rr + atan(fabs(dzx_save/dzy_save)) * 70); 
  1019. //            i = i%(NMAX/2);
  1020. //            if ((int)rr%2 == 1)
  1021. //                i += (NMAX/2);
  1022.             return;
  1023.  
  1024.         case 142:
  1025.             i = (int)(rr + atan(fabs(dzx_save/dzy_save)) * 50); 
  1026.             return;
  1027.  
  1028.         case 143:
  1029.             i = (int)(rr/3 + atan(fabs(dzx_save/dzy_save)) * 30); 
  1030.             if ((int)atan(fabs(dzx_save/dzy_save))%2 == 0 && i < (NMAX/2))
  1031.                 i += (NMAX/2);
  1032.             return;
  1033.  
  1034.         case 144:
  1035.             i = (int) (atan(fabs(z.imag()/z.real())) * 13 + (rr * 10));
  1036.                 
  1037.       return;
  1038.  
  1039.         case 145:
  1040. //            if (rr > 1)    
  1041. //                i = (int) (atan(fabs(z.imag()/z.real())) * nBay100 * 1
  1042. //                    * (rr * nBay1000 * 1));
  1043.       
  1044. //                i = (int) (atan(fabs(z.imag()/z.real())) * nBay100 * 1
  1045. //                    + (rr * nBay1000 * 5));
  1046.                         
  1047. //            if (rr > 1)    
  1048. //                i += (int) (rr * nBay100 * 10 *
  1049. //                    (atan(fabs(z.imag()/z.real())) * nBay100));
  1050. //            else
  1051. //       i = 0;    
  1052.  
  1053.             if (rr >= 1)    
  1054.                 i = (int) (atan(fabs(z.imag()/z.real())) * nBay100 * 5
  1055.                     + (rr * nBay1000 * 5));
  1056.             else
  1057.         i = 0;      
  1058.             
  1059.       return;
  1060.  
  1061.         case 146:
  1062.             i = (int)rr; 
  1063.             return;
  1064.  
  1065.         case 147:
  1066.             i = (int)(rr + atan(fabs(dzx_save/dzy_save)) * (10+nBay1000)); 
  1067.             return;
  1068.  
  1069.     }
  1070.  
  1071.     if (bMFilter)
  1072.     {
  1073.         if (i >= NMAX)
  1074.             i = NMAX - 1;
  1075.         
  1076.         // Set up to Calculate the Fractal Dimension
  1077.         if (i < NMAX)
  1078.         {
  1079.             // Finish out the array with the last value
  1080.             for (i2 = i ; i2<NMAX ; i2++)
  1081.             {
  1082.                 pXTemp[i2] = z.real();
  1083.                 pYTemp[i2] = z.imag();
  1084.             }
  1085.         }
  1086.         
  1087.         // Initialize the mean with zero
  1088.         x_mean = 0;
  1089.         y_mean = 0;
  1090.  
  1091.         // Initialize the range    with the first value
  1092.         x_rmin = pXTemp[0];    // x range min
  1093.         x_rmax = pXTemp[0];    // x range max
  1094.         y_rmin = pYTemp[0];    // y range min
  1095.         y_rmax = pYTemp[0];    // y range max
  1096.  
  1097.         switch(nFilter)
  1098.         {
  1099.             case 29:            // Fractal Dimension
  1100.                 if (i < 3)
  1101.                     return;
  1102.  
  1103.                 if (bQuickMode)
  1104.                 {
  1105.                     iQuickAbort = NMAX;
  1106.                     JMAX = NMAX = i;
  1107.                 }
  1108.                 
  1109.                 for (i2 = 0 ; i2<NMAX ; i2++)
  1110.                 {
  1111.                     // Get sum x and y
  1112.                     x_mean += pXTemp[i2];                                            
  1113.                     y_mean += pYTemp[i2];
  1114.  
  1115.                     // Get min x
  1116.                     if (pXTemp[i2] < x_rmin)
  1117.                         x_rmin = pXTemp[i2];
  1118.  
  1119.                     // Get max x
  1120.                     if (pXTemp[i2] > x_rmax)
  1121.                         x_rmax = pXTemp[i2];
  1122.  
  1123.                     // Get min y
  1124.                     if (pYTemp[i2] < y_rmin)
  1125.                         y_rmin = pYTemp[i2];
  1126.  
  1127.                     // Get max y
  1128.                     if (pYTemp[i2] > y_rmax)
  1129.                         y_rmax = pYTemp[i2];
  1130.                 }
  1131.  
  1132.                 x_mean = x_mean / NMAX;
  1133.                 y_mean = y_mean / NMAX;
  1134.  
  1135.                 x_std = 0;
  1136.                 y_std = 0;
  1137.  
  1138.                 for (i2 = 0 ; i2<NMAX ; i2++)
  1139.                 {
  1140.                     x_std = x_std + (pXTemp[i2] - x_mean)*(pXTemp[i2] - x_mean); 
  1141.                     y_std = y_std + (pYTemp[i2] - y_mean)*(pYTemp[i2] - y_mean); 
  1142.                 }
  1143.  
  1144.                 x_std = sqrt(x_std / (NMAX-1));        // Standard Deviation x
  1145.                 y_std = sqrt(y_std / (NMAX-1));        // Standard Deviation y
  1146.  
  1147.                 // Initialize complex variables
  1148.                 cx_std = cmplx(x_std, 0);    // Complex Standard deviation    x
  1149.                 cy_std = cmplx(y_std, 0);    // Complex Standard deviation    y
  1150.  
  1151.                 cFDx = cmplx (1, 0);                // Complex Fractal Dimension x
  1152.                 cFDy = cmplx (1, 0);                // Complex Fractal Dimension y
  1153.  
  1154.                 cRng_x = cmplx (x_rmax - x_rmin, 0);    // Complex Range x
  1155.                 cRng_y = cmplx (y_rmax - y_rmin, 0);    // Complex Range y
  1156.  
  1157.                 cNMAX = cmplx (NMAX, 0); 
  1158.                 
  1159.                 //////////////////////////////////////////////////////////////
  1160.                 // Fractal Dimension Equation ...
  1161.                 // Ju = Upper Jaenisch coefficient 
  1162.                 // Initialize Ju=1;
  1163.                 // Iterate Ju = [log(R/S*N*J]/[log(1/N)] until abs(Ju+1)-Ju < .0001.
  1164.                 // if ( 0  <= J <= .5,) D(j) = 1/(1-Ju)
  1165.                 // if (.5  <=    J <= 1.0) D(j) = 1/Ju.
  1166.                 ///////////////////////////////////////////////////////////////
  1167.  
  1168.                 //////////////////////////////////////////////////////////////
  1169.                 // Real
  1170.                 ////////////////////////////////////////////////////////////                
  1171.                 dm = 99;
  1172.                 da = 0;
  1173.                 nDIter_x = 0;
  1174.                 denominator = log((double)(1.0/(double)NMAX));
  1175.                 //while (fabs(da - dm) > .0001 && nDIter_x++ < NMAX-2)
  1176.                 while (fabs(da - dm) > .0001 && nDIter_x++ < 10000)
  1177.                 {
  1178.                     // Calculate Fractal Dimension (real)
  1179.                     da = dm;
  1180.                     cFDx = (cRng_x/(cx_std*cNMAX*cFDx)).clog()/denominator;
  1181.                     dm = cFDx.real();
  1182.                     if (nDIter_x == 1)
  1183.                         dFDx_0 = 2.0 - dm;
  1184.                 }
  1185.  
  1186.                 dFDx = cFDx.real();
  1187.  
  1188.                 if (dFDx <= .5)
  1189.                     dFDx = 1 / (1 - dFDx);
  1190.                 else
  1191.                     dFDx = 1 / dFDx;
  1192.  
  1193.                 ////////////////////////////////////////////////////////////                
  1194.                 // Imaginary
  1195.                 ////////////////////////////////////////////////////////////                
  1196.                 dm = 99;
  1197.                 da = 0;
  1198.                 nDIter_y = 0;
  1199.                 denominator = log((double)(1.0/(double)NMAX));
  1200.                 //while (fabs(da - dm) > .0001 && nDIter_y++ < NMAX-2)
  1201.                 while (fabs(da - dm) > .0001 && nDIter_y++ < 10000)
  1202.                 {
  1203.                     // Calculate Fractal Dimension (real)
  1204.                     da = dm;
  1205.                     cFDy = (cRng_y/(cy_std*cNMAX*cFDy)).clog()/denominator;
  1206.                     dm = cFDy.real();
  1207.                     if (nDIter_y == 1)
  1208.                         dFDy_0 = 2.0 - dm;
  1209.                 }
  1210.  
  1211.                 dFDy = cFDy.real();
  1212.  
  1213.                 if (dFDy <= .5)
  1214.                     dFDy = 1 / (1 - dFDy);
  1215.                 else
  1216.                     dFDy = 1 / dFDy;
  1217.                 ////////////////////////////////////////////////////////////                
  1218.                 
  1219.                 ////////////////////////////////////////////////////////////
  1220.                 // Apply Fractal Dimension to a color
  1221.                 
  1222.                 switch (nFDOption)
  1223.                 {
  1224. //                    case 0:
  1225. //                        // Addition
  1226. //                        jrw = (int) ((dFDx+dFDy) * 100 * nBay100);
  1227. //                        break;
  1228.                 
  1229.                     case 1:
  1230.                         // Fractal Dimension Real Initial (FD_0)
  1231.                         jrw = (int) (dFDx_0 * 100.0 * nBay100);
  1232.                         break;
  1233.                 
  1234.                     case 2:
  1235.                         // Fractal Dimension Imginary    Initial (FD_0)
  1236.                         jrw = (int) (dFDy_0 * 100.0 * nBay100);
  1237.                         break;
  1238.                 
  1239.                     case 3:
  1240.                         // (FD_0) Addition
  1241.                         jrw = (int) ((dFDx_0+dFDy_0) * 100 * nBay100);
  1242.                         break;
  1243.  
  1244.                     case 4:    
  1245.                         // (FD_0) Sum of Squares
  1246.                         jrw = (int) ((dFDx_0*dFDx_0+dFDy_0*dFDy_0) * 100 * nBay100);
  1247.                         break;
  1248.  
  1249.                     case 5:
  1250.                         // (FD_0) Multiplication
  1251.                         jrw = (int) (dFDx_0 * dFDy_0 * 100 * nBay100);
  1252.                         break;
  1253.                     
  1254.                     case 6:
  1255.                         // Fractal Dimension Real
  1256.                         jrw = (int) (dFDx * 100.0 * nBay100);
  1257.                         break;
  1258.                 
  1259.                     case 7:
  1260.                         // Fractal Dimension Imginary
  1261.                         jrw = (int) (dFDy * 100.0 * nBay100);
  1262.                         break;
  1263.                 
  1264.                     case 8:
  1265.                         // FD Addition
  1266.                         jrw = (int) ((dFDx+dFDy) * 100 * nBay100);
  1267.                         break;
  1268.  
  1269.                     case 9:    
  1270.                         // FD Sum of Squares
  1271.                         jrw = (int) ((dFDx*dFDx+dFDy*dFDy) * 100 * nBay100);
  1272.                         break;
  1273.  
  1274.                     case 10:
  1275.                         // FD Multiplication
  1276.                         jrw = (int) (dFDx * dFDy * 100 * nBay100);
  1277.                         break;
  1278.                                     
  1279.                     default:
  1280.                         AfxMessageBox("FD Error");
  1281.                         break;
  1282.                 }                
  1283.  
  1284.                 if (bQuickMode)
  1285.                 {
  1286.                     //i = NMAX + abs(jrw);
  1287.                     i = abs(jrw);
  1288.                     NMAX = iQuickAbort;
  1289.                     if (i >= NMAX) i %= (NMAX-1);
  1290.                 }
  1291.                 else
  1292.                 {
  1293.                     i = abs(jrw);
  1294.                     if (i >= NMAX) i %= (NMAX-1);
  1295.                 }
  1296.  
  1297.                 return;
  1298.  
  1299.                 break;
  1300.  
  1301.             default:            // Standard deviation by default
  1302.                 for (i2 = 0 ; i2<i ; i2++)
  1303.                 {
  1304.                     // Get sum x and y
  1305.                     x_mean += pXTemp[i2];                                            
  1306.                     y_mean += pYTemp[i2];
  1307.  
  1308.                     // Get min x
  1309.                     if (pXTemp[i2] < x_rmin)
  1310.                         x_rmin = pXTemp[i2];
  1311.  
  1312.                     // Get max x
  1313.                     if (pXTemp[i2] > x_rmax)
  1314.                         x_rmax = pXTemp[i2];
  1315.  
  1316.                     // Get min y
  1317.                     if (pYTemp[i2] < y_rmin)
  1318.                         y_rmin = pYTemp[i2];
  1319.  
  1320.                     // Get max y
  1321.                     if (pYTemp[i2] > y_rmax)
  1322.                         y_rmax = pYTemp[i2];
  1323.                 }
  1324.  
  1325.                 if (i != 0)
  1326.                 {
  1327.                     x_mean = x_mean / (double)i;
  1328.                     y_mean = y_mean / (double)i;
  1329.                 }
  1330.  
  1331.                 x_std = 0;
  1332.                 y_std = 0;
  1333.  
  1334.                 for (i2 = 0 ; i2<i ; i2++)
  1335.                 {
  1336.                     x_std = x_std + (pXTemp[i2] - x_mean)*(pXTemp[i2] - x_mean); 
  1337.                     y_std = y_std + (pYTemp[i2] - y_mean)*(pYTemp[i2] - y_mean); 
  1338.                 }
  1339.                                 
  1340.                 if (i>1)
  1341.                 {
  1342.                     x_std = (x_std / (i-1));
  1343.                     y_std = (y_std / (i-1));
  1344.                 }
  1345.                 else
  1346.                 {
  1347.                     x_std = 0;
  1348.                     y_std = 0;
  1349.                 }
  1350.  
  1351.                 if (x_std == 0 && y_std == 0)
  1352.                     jrw = 0;
  1353.                 else
  1354.                     jrw = (int) (log((x_std + y_std)/2) * 10 * nBay100);
  1355.  
  1356.                 i = abs(jrw);
  1357.                                         
  1358.                 if (i >= NMAX) i %= (NMAX-1);
  1359.                 return;
  1360.         }
  1361.     }
  1362.  
  1363.     if (i<NMAX)
  1364.         i += jrw;
  1365.     else
  1366.         i = jrw;
  1367.  
  1368.     i = abs(i);
  1369.     if (i >= NMAX) i %= (NMAX-1);
  1370. }     
  1371.